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Abstract 

A new Monte Carlo move for polymer simulations 
is presented. The "wormhole" move is build out 
of rcptation steps and allows a polymer to reptate 
through a hole in space; it is able to completely 
displace a polymer in time N"^ (with N the poly- 
mer length) even at high density. This move can 
be used in a similar way as configurational bias, 
in particular it allows grand canonical moves, it is 
applicable to copolymers and can be extended to 
branched polymers. The main advantage is speed 
since it is exponentially faster in N than configu- 
rational bias, but is also easier to program. 



as CB but it is much faster since it is quadratic in N 
instead of exponential. Based on reptation moves, 
our algorithm is relatively easy to program. Never- 
theless our method is quite different from standard 
reptation since it can do long-range displacements 
and grand-canonical moves, it is also much more 
flexible and can be applied to hetero and branched 
polymers. 

In the next section, we precisely describe this new 
method and prove its correctness, then we show 
with numerical benchmarks that our algorithm is 
much faster than CB and finally we introduce a 
set of possible modifications (in particular grand 
canonical moves). 



1 Introduction 

Polymer systems are very common (DNA, proteins, 
plastics, . . . ) and theoretical studies are quite 
difficult. That is why efficient numerical simula- 
tion techniques are required, either to test theo- 
retical approximations or to compare models with 
experiments.^ Unfortunately polymer simulations 
are very slow because they involve big time and 
length scales. A lot of smart Monte Carlo moves 
have been devised to reduce the relaxation time, 
of such systems, in particular the pivot algorithmcl 
for isolated polymer chains, the rjeptation method 
(or "slithering snake" algorithm) J3 ccmfigurational 
bias (CB)qI3 and extensions around it.ErQ Moreover 
general Monte Carlo methods such as exchange 
Monte Carlo, (or parallel tempering) histogram 
re-weightingJ or "go with the winner" simulationsli2l 
are also used. 

Here we propose a new kind of Monte Carlo move 
for polymer simulations. This move, called "worm- 
hole" (WH) move, has more or less the same effect 



2 The Algorithm 

We first roughly describe our algorithm. Wc con- 
sider a linear polymer of length iV in a given 
fixed environment (temperature, other polymers, 
solvent, . . . ). The WH Monte Carlo move opens 
a hole in space (the wormhole) and tries to make 
the polymer reptate through it (see figure |]) . The 
move is composed of a lot of elementary reptation 
steps which are individually rejected or accepted 
according to standard Metropolis. This goes on 
as long as the polymer is not in one piece again. 
Each of the reptation step is randomly backward 
or forward and the polymer thus does a discrete 
one-dimensional random walk. Three different sce- 
narios may happen: 

1. The very first step fails, the hole is not opened 
and nothing happens. We call it a failed move. 

2. The polymer reptates through the hole as 
sketched in figure |^; the random walk has cov- 
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Figure 1: A complete WH move. Note how the 
monomers move as they come and go through the 
hole. 



ered a distance N. The polymer is in a com- 
pletely new location and new configuration. 
This is a complete move. 

3. The polymer enters the hole but comes out on 
the same side; the random walk got back to 
its starting point. The monomers which got 
through the hole and came back have moved 
and the others are unchanged. This is a partial 
move. 

Let us now enter into the details and describe 
exactly what the algorithm does. WH proceeds as 
follow: 

1. Wormhole drilling step: Randomly choose one 
end-monomer (with probability 1 /2) and try to 
move it to a uniformly random position, with 
the old bond broken and a virtual one created 
to the other end of the polymer (the dashed 
line in figure |l|). Proceed to step 3. 

2. Standard reptation step: Randomly choose 
one end-monomer and try to move it to the 
other end of the polymer connecting it with 
a uniformly random bond. For this move, we 
do as if the virtual bond was a normal one, so 
that the polymer has only two ends. 



3. End test: If the polymer is in two pieces, pro- 
ceed to step 2. Otherwise the WH move is 
finished and is accepted with probability one. 

At step 1 and 2, the move is accepted with the 
Metropolis probability: 



Fa/(A£;) =min(l, 



(1) 



where /3 = l/ksT and AE is the energy difference 
generated by the trial move. Note that the virtual 
bond may be given a constant energy to improve 
the acceptance rate of the first step. 

The claim is that WH is a correct Monte Carlo 
move: it overall respects the detailed balance con- 
dition: 

P{C)P{C -> C") = P{C')P{C' C), (2) 

where P{C C) is the transition probability from 
configuration C to C" and P(C) = e-'^^C^) /Z is the 
Boltzmann weight. It is easy to check that the new 
configuration is a valid one : (i) there is neither 
hole nor virtual bond any longer, (ii) the nature of 
the monomers and their order along the chain are 
the same as in the old configuration (that is why 
it works for hetero-polymers, in contrast to simple 
reptation). 

If the first step is rejected, equation |^ is trivially 
fulfilled. Otherwise, n elementary steps are per- 
formed and the system goes through a sequence of 
configurations Cq, • • • , C„ with energies Eq, - ■ ■ , En- 
Each elementary step Si = Ci-i — > Ci takes the 
system from one configuration to the next (with 
Ci-i = Ci if the step is rejected). The probability 
of the sequence S = {Si, • • • , Sn) can be written 

n 

P{S)^l[P{S.). (3) 

1=1 

To each sequence corresponds an opposite sequence 
S = {Sm ■ ■ ■ , Si), made of the opposite steps Si = 
Ci — *■ Ci-i, which takes the system from C„ to Co 
and we have 



P{S) 



n 



For all accepted intermediate steps (1 < i < n), we 
have: 



P{S^) 



P{S) f \ P{S.) 

crmediate steps 

\PB{S^)PM{E^- E,-i) 

\Pb{Si)Pm{E,-i - E,) 



(4) 



(5) 
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where Pb is the uniform probability of choosing 
a bond and the factor 1/2 comes from the choice 
of the reptation direction. The rejected steps fulfil 
Si = Si and cancel. They thus also fulfil equation||. 
On the other hand, for i ~ 1, 



PjSi) ^ Pp{S^) ^ 
P(^) Pb(^)' 



-f3{Ei-Eo) 



and for i 



P{Sn) _ PsiSn) -(}{E^-E. 

p(s::) Pp(5;) 



(6) 



(7) 



with Pp the uniform probability of choosing a po- 
sition within the simulation box. Hence we get 



PjS) 
PiS) 



-P(E„~Eq) 



(8) 



Since P{C C) is the sum of ^(5*) over all pos- 
sible sequence S with Co = C and Cn = C , we 
get equation ^. This proves the correctness of our 
algorithm. 

Finally, concerning ergodicity, it is clear that this 
algorithm is able to move a polymer to any place 
where the required space exists beforehand. But, 
as for reptation and CB, one cannot for example 
ensure global ergodicity on a very densely occupied 
lattice. Nevertheless, as shown in the next section, 
WH appears to perform quite well even at relatively 
high density. 

3 Efficiency 

Now that the algorithm is proven correct, let us 
see how good it is. Usually the main difficulty for 
a Monte Carlo algorithm for polymer is to relax 
the configuration of the polymers, which is more 
an entropy problem than an energy problem. That 
is why we test the speed of this neaz. algorithm on 
the bond fluctuation model (BFM)tiHl3 at infinite 
temperature. The BFM is a lattice polymer model 
where each monomer occupies 8 sites (2 x 2 x 2) of 
a cubic lattice and monomers cannot overlap (hard 
core repulsion). The bond vectors are restricted to 
a set of 108 possibilities with lengths 2, -s/S, \/5, 3 
and VTO. The interaction range is set to a/6 which 
means that only monomers in direct contact inter- 
act. 




Figure 2: Comparison of the conformations pro- 
duced by WH and by standard reptation for an 
isolated polymer of size iV = 50 at T = 2. Both 
figures show two curves superimposed, but they 
are indistinguishable. Top: the integrated distribu- 
tions of cos 6 with 9 the angle between two succes- 
sive bonds. Bottom: the integrated distributions 
of sin (p with (p the (torsion) angle between a bond 
and the plane defined by the two preceeding bonds. 
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Figure 3: The number of accepted elementary steps 
required to achieve a complete WH move is 0{N'^). 



To check the correctness of the implementation, 
we used our program to simulate isolated chains of 
different lengths with attractive interactions near 
the 8-temperature (namely T = 2). This is the 
most sensitive point to check that the conforma- 
tion of the chain is correctly sampled, since it is 
where the fluctuations and the temperature depen- 
dence are the largest. First we reproduced the mea- 
surern£nts of the 0-temperature done by Wilding 
et alJEi and found fully compatible results (namely 
Te = 2.03(3) to be compared with Te = 2.02(2)). 
Moreover we compared the distributions of the 
bond-bond and torsion angles with the ones given 
by standard reptation. The results are shown on 
figure ||. Both algorithms give the same results 
(within statistical error bars, here roughly 0.1%) 
and the curves, on top of one another, are indistin- 
guishable. Note that the roughness of the curves is 
no noise, it is due to the discreteness of the model. 

All the following benchmarks are made using a 
mixture of polymers of length N at density 1/2 in 
the BFM at infinite temperature. For these param- 
eters, our computer (a Pentium III at 500 MHz) 
achieve 1.3 10^ elementary steps per second with 
an acceptance rate r ~ 8.5%. Note that our imple- 
mentation also works at finite temperature and is 
thus slower than an implementation optimised for 
the infinite temperature case. 

First, we have claimed that our algorithm is able 
to completely displace a polymer in time N'^. This 
is a reasonable assumption since it is the time re- 
quired by an unbiased random walk to cover a dis- 



Figure 4: CPU time ratios between CB and WH as 
a function of N. Crosses are time ratios for com- 
pletely displacing a polymer (time to accept a CB 
move/time to do a complete WH move). Diamonds 
are relaxation time ratios (see equation |o|) . 

tance N and this corresponds to a complete WH 
move. But this must be tested numerically since 
the random walk actually done is strongly infiu- 
enced by density fluctuations. We note m the mean 
number of elementary steps required to achieve one 
complete WH move (many failed and partial moves 
may also happen in between) and mr is the cor- 
responding number of accepted elementary steps. 
The values of mr/N'^ are shown as a function of N 
in flgure ^, they clearly show that for large N, we 
have 

m = -iV2, (9) 
r 

with a ~ 800. To achieve a complete move on our 
computer, we thus need around 7. 10~^A^^ seconds 
for a long polymer (and less for a smaller one). As 
could be expected, the value of a decreases at lower 
density, for example at density 0.1, a ~ 1.8 (and 
r ~ 67.5%). 

What our algorithm can do is very similar to CB, 
so it is interesting to do the comparison. The CB 
with which we compare tries to erase a polymer 
completely and to recreate it somewhere else, for 
each monomer it checks the 108 possible bonds and 
chooses an allowed one (not already occupied) and 
computes the associated weight. On our computer, 
we can do 7.10^^ such elementary steps per second 
(one needs N of them for a whole CB move). As 
for WH, our CB implementation also works at fl- 
nite temperature and is thus slower than it should 
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for the infinite temperature case. Note that one 
could also use a CB implementation which does not 
check all the 108 bonds (say 10 or 20). In this case 
the elementary steps are accordingly faster but the 
overall acceptance rate decreases much faster with 
N. 

In figure 0, we show the ratios of time required 
by both algorithm to completely displace a polymer 
(an accepted CB move and a complete WH move). 
From this point of view, the algorithms are more or 
less equivalent at small N, and WH becomes expo- 
nentially better at large N. In this first comparison 
we did not take into account the fact that WH also 
relaxes the system when it does a partial move. So 
we measured the relaxation time of the end-to-end 
vector R(i), namely: 



r = -i=-y' (R(0).R(0)dt 



(10) 



The ratios of those times are also shown on figure |j, 
clearly WH is much faster than CB. On the other 
hand, at a lower density, CB is more efficient for the 
small values of iV (for example at density 0.1 and at 
N = 40, CB is only 3 times slower than WH instead 
of 100 times slower), but CB stays exponentially 
slower at large N. See below how one could take 
advantage of a faster CB for small N. 

Finally this comparison was done on a lattice, 
this gives a big advantage to CB since off-lattice 
CB implementations require a large number of en- 
ergy computations which are much longer than on 
a lattice. Moreover, reptation is hampered by the 
lattice since it can get blocked much more easily 
than off-lattice. We thus think that WH should be 
even better for continuous models. 

Our WH method is reminiscent of the extended 
ensemhlcLjnethod introduced by Escobedo and de 
Pablo.oEj But in their method, the elementary 
steps are CB moves that move k monomers at a 
time. Moreover the intermediate configurations are 
also sampled giving rise to extended ensembles with 
configurations where a polymer is incomplete or 
split into two pieces. Since they relax the whole 
system between each elementary step the time to 
deal with one polymer is 



a'V 



m 



(11) 



with V the number of particles in the system and a' 
and r' the equivalent of a and r in equation ^. Here 



r' oc exp(— fc/fco) and there is an optimal value for 
k (namely k = 2ko). With our parameters /cq — 4.5 
and r' ~ 3%, so even if a' < a, V/k"^ is very large. 

Finally, WH has one feature which may be un- 
desirable. The execution time for one move has 
huge fluctuations (that explains the fluctuations in 
figure A move can be very fast (for example 
a Jailed move) but it can also last a few seconds 
(and exceptionally a few minutes for long chains). 
In a simple simulation that is no problem, since 
only the average time is relevant. But in a parallel 
computation, where each processor waits after the 
slowest one, this can become a real difficulty. A 
first solution is to do a lot of WH moves between 
two synchronisations, thus averaging the fluctua- 
tions. Another one is to set a limit to the number 
of elementary steps which can be done in one move 
(say 10m). If this threshold is exceeded, the move 
is rejected and the polymer is restored to its initial 
state. 



4 Variations 

Let us first see algorithmic refinements to WH. 
First, in the case of an off-lattice simulation one 
may wish to draw the polymer bonds out of a non- 
uniform Pb distribution. In particular one may 
want to draw the bond length around a certain 
favourable value. In such a case it is necessary to 
add a correcting weight to the Metropolis rule to 
recover equation |^ (as done for normal reptation) . 
Then one can check that these correcting weights 
raise the correct overall correction term and equa- 
tion H still holds so the algorithm is still correct. 
More generally, as soon as one wants to modify the 
WH algorithm, one should check that the cancella- 
tions still occur, in particular for i ^ 1 and i = n. 

In the original reptation methocH the direction 
of the reptation is changed only when a move is 
rejected and not drawn randomly at each step as 
we propose here. Nothing prevents to do the same 
in WH (the proof just gets a little more compli- 
cated). At very low density this can improve the 
speed by a factor 2 or more because successive steps 
are done in the same direction. But at high density, 
it changes essentially nothing. 

As explained in the previous section, CB can be 
more efficient at small TV than WH (at low density 
or when we do not check all the 108 bonds at each 
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step). Then, as done in the "extended ensemble" 
method, one can use CB as the elementary step, 
moving k monomers at a time and we then have 
(with notation of the previous section) 



a 

r' \ k 



N 



(12) 



and the CPU time is proportional to km (the opti- 
mal value of k is now k — ko). This is faster than 
the original WH if fcr'/r > K with K the CPU time 
ratio to do one elementary move with both meth- 
ods (here K ~ 19). Note that we have supposed 
that a does not change which is not clear. Finally 
one of the interesting point in WH is the ease of 
programming, this is lost in this variation. 

Let us now see how WH can be extended. First, 
how to do grand canonical moves. These moves can 
be of three different kinds: (i) to move a polymer 
from one simulation box to another one; (ii) to re- 
move one polymer from a box; (iii) to insert a new 
polymer in a box. The first one is very easy to do, 
WH applies directly and one can eventually add a 
chemical potential for each monomer depending on 
the box. For (ii) and (iii) , one proceeds more or less 
as for CB. The algorithm has to be modified, the 
reptation steps are made between the simulation 
box and some kind of black box without geometry 
where the monomer have a given chemical poten- 
tial. The Pb and Pp no longer completely cancel 
(the steps toward the black box do not require a 
random draw) and for an inserting move and the 
corresponding removing move, equation ^ becomes 



^ ^^PiE,,-Eo)p pN- 

PiS) ^ 



(13) 



The new factor is a constant which can be absorbed 
in the chemical potential (exactly as the weight 
with CB). This factor is exactly the probability for 
choosing a particular configuration of the polymer. 
We did not perform tests for the grand canonical 
moves, but we expect them to behave as described 
in the previous section. 

Very often CB is not used to remove and insert a 
whole polymer, but is applied only on a part of it. 
Such a thing can also be done with WH, provided 
we perform a little modification (see figure ^) . At 
the first step of the algorithm, instead of choos- 
ing a random position, the monomer is randomly 
bound to another non-end monomer of the same 








Figure 5: A complete modified WH move. The 
letters represent different kind of monomers in a 
hetero-polymer. Note how in the intermediate 
configurations, monomer F has disappeared and 
monomer C is present twice. 




Figure 6: The order in which the monomers enter 
the hole (left), and leave it (right) for a branched 
polymer. 



polymer. After this, WH proceeds as before until 
the polymer is linear again. Note that in the case of 
a hetero-polymer, the monomer that goes through 
the hole may be replaced by another one, so that 
the polymer has kept its structure at the end (see 
figure I). This modified move is probably less effi- 
cient than the original one for the same number of 
monomers, since the new growing branch is going to 
encounter an already densely occupied region (the 
old branch being also there). Moreover the origi- 
nal WH is also able to move parts of the polymer 
(with partial moves), so it is perhaps better to use 
the original move. The modified move can in any 
case be useful for branched polymers. 

Finally one can also extend the original WH to 
branched polymers. The basic idea stays the same: 
a random walk through a set of non-valid configu- 
rations where the polymer is split into two pieces. 
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But it cannot use simple reptation because of the 
branching points. One way to achieve this is to 
make the monomers get out of the hole in a dif- 
ferent order as they enter it. See figure ^ for an 
example. As in the previous paragraph, during 
the move, certain monomers may be present twice 
whereas others are no longer there. 

5 Discussion and Conclusions 

The wormhole move (WH) presented here is a new 
kind of Monte Carlo move for polymer simulations. 
It can completely displace a polymer in time 0{N'^) 
even at high density (1/2 in the bond fluctuation 
model). It is based on reptation steps and is thus 
quite easy to program. It is nevertheless very flexi- 
ble since it works on hetero-polymers, allows grand 
canonical moves (insertions and deletions of poly- 
mers) and can be adapted to branched polymers. 

What is this new method good for ? Essentially, 
it is meant to replace configurational bias (CB). 
The wormhole move can in fact do whatever CB 
does, and it does it faster. As CB, it can destroy, 
create or displace a polymer, but CB works in time 
0{e^^) whereas WH works in time 0{N'^). Thus, 
exception made of the case of short polymers at 
low density, it is always better to use WH instead 
of CB. In the case of short polymers at low density, 
it is probably better to use WH as well, since this 
case is anyway an easy one which requires short 
computation time and WH is relatively easier to 
program. 

On the other hand, the aim of WH is not to re- 
place the standard reptation method or to com- 
pete with it. It simply does something different. 
Typically standard reptation is to be preferred in 
the case of linear homo-polymers with small den- 
sity fluctuations (i.e. when it is not important to 
be able to rapidly displace a polymer between two 
distant positions). In the other cases, where one 
would have used CB, one should use WH. There 
are essentially two main cases: (i) if reptation is 
not applicable, in particular to do grand canonical 
moves or to simulate hetero-polymers; (ii) in the 
case of large density fluctuations (typically when 
more than one phase is present) where one wants 
to rapidly relax the density by displacing polymers 
over large distances. 

We expect this new method to become widely 



used since it is easy to program and very efficient. 
At the moment, it is ojsed to simulate dense random 
copolymer mixtures .113 

The author acknowledges M. Miiller, L. G. Mac- 
Dowell, A. Yethiraj and K. Binder for fruitful dis- 
cussions. 
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